% This file simulates the path of household inflation expectation
% given a monetary policy that increases the inflation target by 1%
% Simulations are done using model-implied parameters

% hh.psi = psi_{t|T}
% t1, t2 = time of interest
% hor = length of IRF

clear;
close all;
paths;

load('soln.mat');
load('data.mat');

Nhh = 2;
hor = 20;

data.alpha = alpha(1,:)'; data.alpha_nc = alpha(3,:)'; data.alpha_cd = alpha(4,:)';
data.phi = phi0(1,:)'; data.phi_nc = phi0(3,:)'; data.phi_cd = phi0(4,:)';
data.px1_actual = pi_actual'; data.px1_nc = infl_e(2,:)'; data.px1_cd = infl_e(3,:)';

hh.K = NaN(3,1,hor,Nhh); hh.P = NaN(3,3,hor+1,Nhh);
hh.psi = NaN(3,hor+1,Nhh);
cc = 0.1;
simul = 1000; exp_nc = NaN(hor+1, simul);

irf_inflation = NaN(Nhh+1, hor+1, 2);
irf_expectation = NaN(hor+1, 1, Nhh+1);

event = 1;

for period = [20122, 19792]
    
    ttt = find(data.YYYYQ == period);

    % MP shock
    PI = data.px1_nc(ttt);
    ALPHA = data.alpha(ttt);
    PHI = 1/(1+cc*exp(-data.phi(ttt)));
    CC = ALPHA/(1-PHI);


    if period <= 19900
        shock = (PI-1)*(1-PHI);
    else
        shock = (PI + 1)*(1-PHI);
    end
    
    
    s_a = NaN(Nhh+1, hor+1); 
    s_a(:,1) = shock +  PHI * PI;
    
    for t = 1:hor
        s_a(:,t+1) = shock + PHI * s_a(:,t);  % inflation
    end
    
    irf_inflation(:,:,event) = s_a;
    
    s_b = repmat(s_a(2,:)', [1, simul]) + normrnd(0, soln.hh_sd_w(3), [hor+1, simul]); % signal to non-college
    
    
    % rational

    d_PHI = cc*exp(data.phi(ttt))/(cc+exp(data.phi(ttt)))^2;
    [~,spf1_P] = hh_K(soln.sd_eta^2, soln.sd_nu1^2, soln.sd_nu2^2, 0, diag([1,1,1]), PHI, PI, d_PHI);
    for iii = 1:20
        [~,spf1_P] = hh_K(soln.sd_eta^2, soln.sd_nu1^2, soln.sd_nu2^2, 0, spf1_P, PHI, PI, d_PHI);
    end
        
    spf1.P(:,:,1) = spf1_P;
    spf1.psi(:,1) = [PI, ALPHA, data.phi(ttt)];
    
    for t = 1:hor
            
        phi = 1/(1+cc*exp(-spf1.psi(3,t)));
        d_phi = cc*exp(spf1.psi(3,t))/(cc+exp(spf1.psi(3,t)))^2;
        [spf1.K(:,:,t) , spf1.P(:,:,t+1)] = hh_K(soln.sd_eta^2, soln.sd_nu1^2, soln.sd_nu2^2, 0, spf1.P(:,:,t), phi, spf1.psi(1,t), d_phi);
        psi_t = spf1.psi(:,t) + spf1.K(:,:,t) * (s_a(1,t) - spf1.psi(1,t));
        spf1.psi(1,t+1) = psi_t(2) + psi_t(1)/(1+cc*exp(-psi_t(3)));
        spf1.psi(2:3,t+1) = psi_t(2:3);
    
    end
    irf_expectation(:,event,Nhh+1) = spf1.psi(1,1:t+1)';
    
    % college
    n = 2;
    hh.P(:,:,1,n) = est.hh_P(:,:,ttt,n+1);
    hh.psi(:,1,n) = [data.px1_cd(ttt); data.alpha_cd(ttt); data.phi_cd(ttt)];
    
    for t = 1:hor
            
        phi = 1/(1+cc*exp(-hh.psi(3,t,n)));
        d_phi = cc*exp(hh.psi(3,t,n))/(cc+exp(hh.psi(3,t,n)))^2;
        [hh.K(:,:,t,n) , hh.P(:,:,t+1,n)] = hh_K(soln.sd_eta^2, soln.sd_nu1^2, soln.sd_nu2^2, soln.hh_sd_w(n+1)^2, hh.P(:,:,t,n), phi, hh.psi(1,t,n), d_phi);
        psi_t = hh.psi(:,t,n) + hh.K(:,:,t,n) * (s_a(n+1,t) - hh.psi(1,t,n));
        hh.psi(1,t+1,n) = psi_t(2) + psi_t(1)/(1+cc*exp(-psi_t(3)));
        hh.psi(2:3,t+1,n) = psi_t(2:3);
    
    end
    
    irf_expectation(:,event,n) = hh.psi(1,:,n);
    
    % non-college
    n = 1;
    hh.P(:,:,1,n) = est.hh_P(:,:,ttt,n+1);
    hh.psi(:,1,n) = [data.px1_nc(ttt); data.alpha_nc(ttt); data.phi_nc(ttt)];
    
    for count = 1:simul
        
        for t = 1:hor
            phi = 1/(1+cc*exp(-hh.psi(3,t,n)));
            d_phi = cc*exp(hh.psi(3,t,n))/(cc+exp(hh.psi(3,t,n)))^2;
            [hh.K(:,:,t,n) , hh.P(:,:,t+1,n)] = hh_K(soln.sd_eta^2, soln.sd_nu1^2, soln.sd_nu2^2, soln.hh_sd_w(n+1)^2, hh.P(:,:,t,n), phi, hh.psi(1,t,n), d_phi);
            psi_t = hh.psi(:,t,n) + hh.K(:,:,t,n) * (s_b(t,count) - hh.psi(1,t,n));
            hh.psi(1,t+1,n) = psi_t(2) + psi_t(1)/(1+cc*exp(-psi_t(3)));
            hh.psi(2:3,t+1,n) = psi_t(2:3);
        end
        exp_nc(:,count) = hh.psi(1,:,1);
    
    end
    irf_expectation(:,event,n) = median(exp_nc,2);
    


    xline = 0:1:8;
    figure;
    plot(xline, irf_expectation(1:9,1,1), ':', 'LineWidth',3); hold on;
    plot(xline, irf_expectation(1:9,1,2), '-.', 'LineWidth',3); 
    plot(xline, irf_inflation(1,1:9), '-', 'LineWidth',3, 'Color', [0.4940 0.1840 0.5560]);
    legend('boxoff');
    legend({'non-college expected inflation', 'college expected inflation', 'inflation target'},'fontsize', 14, 'Location','northeast');box off;
    ylabel('Inflation', 'fontsize', 16); 
    xlabel('Quarters', 'fontsize', 14); 
    fontsize(gcf,scale=1.5);
    style = hgexport('factorystyle'); style.Color = 'gray';
    figname = sprintf('figures/mp_%.0f.eps', period);
    hgexport(gcf, figname, style);

end




